module my_toolbox
    
    contains
    
    ! ####################################################################################################################
    ! mean of a vector ###################################################################################################
    ! ####################################################################################################################
    subroutine mean(vec,output)
        
        implicit none
        
        double precision, intent(in)  :: vec(:)
        double precision, intent(out) :: output
        
        output = sum(vec) / size(vec)
        
    end subroutine mean
        
    ! ####################################################################################################################
    ! mean of a matrix ###################################################################################################
    ! ####################################################################################################################
    subroutine mean2(matrix,output)
    
        implicit none
        
        double precision, intent(in)  :: matrix(:,:)
        double precision, intent(out) :: output
        double precision :: total_sum
 
        total_sum = sum(sum(matrix,1))
        output = total_sum / ( size(matrix,1) * size(matrix,2) )
        
    end subroutine mean2
        
    ! ####################################################################################################################
    ! mean of a 3-dimensional matrix #####################################################################################
    ! ####################################################################################################################
    subroutine mean3(matrix,output)
        
        implicit none
        
        double precision, intent(in)  :: matrix(:,:,:)
        double precision, intent(out) :: output
        double precision :: total_sum
 
        total_sum = sum(sum(sum(matrix,3),2))
        output = total_sum / ( size(matrix,1) * size(matrix,2) * size(matrix,3) )
        
    end subroutine mean3
        
    ! ####################################################################################################################
    ! mean, variance, standard deviation of a series #####################################################################
    ! ####################################################################################################################
    subroutine std_dev(vec,mean,var,std)
        
        implicit none
        
        double precision, intent(in)  :: vec(:)
        double precision, intent(out) :: mean
        double precision, intent(out) :: var
        double precision, intent(out) :: std
        
        mean = sum(vec) / size(vec)
        var  = sum( (vec - mean)**2d0 ) / size(vec)
        std  = sqrt(var)
        
    end subroutine std_dev

    ! ####################################################################################################################
    ! linear interpolation ###############################################################################################
    ! ####################################################################################################################
    subroutine linear_interp(x, xi, val, interp)    ! Note: works only for evenly-spaced grid!

        implicit none

        double precision, intent(in)    :: x(:)     ! array of grid points
        double precision, intent(in)    :: xi       ! value on the grid line
        double precision, intent(in)    :: val(:)   ! array of values corresponding to grid points
        double precision, intent(out)   :: interp   ! output, interpolated value

        double precision :: step                    ! step size of grid points
        double precision :: x_gap                   ! gap between xi and the lower bound of x
        integer :: loc                              ! index for largest grid point among less than x
        double precision :: ratio                   ! ratio of xi between x(loc) and x(loc+1)

        integer nx
        nx = size(x)                                ! number of grid points

        step   = x(2) - x(1)                                     
        x_gap  = xi - x(1)                                       
        loc    = min( nx - 1, max(1, floor(x_gap/step) + 1) )    
        ratio  = (xi - x(loc)) / (x(loc+1) - x(loc))             
        interp = val(loc) + ratio * ( val(loc+1) - val(loc) )

    end subroutine linear_interp

    ! ####################################################################################################################
    ! 2-dimension linear interpolation ###################################################################################
    ! ####################################################################################################################
    subroutine linear_interp2(x, xi, y, yi, val, interp)   ! Note: works only for evenly-spaced grid!

        implicit none

        double precision, intent(in)    :: x(:)     ! array of grid points
        double precision, intent(in)    :: xi       ! value on the grid line
        double precision, intent(in)    :: y(:)     ! array of grid points
        double precision, intent(in)    :: yi       ! value on the grid line
        double precision, intent(in)    :: val(:,:) ! array of values corresponding to grid points (x,y)
        double precision, intent(out)   :: interp   ! output, interpolated value
        
        double precision :: u1,u2                   ! intermediate values in interpolation
        double precision :: step_x                  ! step size of grid points
        double precision :: gap_x                   ! gap between xi and the lower bound of x
        integer :: loc_x                            ! index for largest grid point among less than x
        double precision :: ratio_x                 ! ratio of xi between x(loc) and x(loc+1)
        double precision :: step_y                  ! step size of grid points
        double precision :: gap_y                   ! gap between xi and the lower bound of x
        integer :: loc_y                            ! index for largest grid point among less than x
        double precision :: ratio_y                 ! ratio of xi between x(loc) and x(loc+1)

        integer nx, ny
        nx = size(x)                                ! number of grid points
        ny = size(y)                                ! number of grid points

        step_x  = x(2) - x(1)                                     
        gap_x   = xi - x(1)                                       
        loc_x   = min( nx - 1, max(1, floor(gap_x/step_x) + 1) )    
        ratio_x = (xi - x(loc_x)) / (x(loc_x+1) - x(loc_x))             

        step_y  = y(2) - y(1)                                     
        gap_y   = yi - y(1)                                       
        loc_y   = min( ny - 1, max(1, floor(gap_y/step_y) + 1) )    
        ratio_y = (yi - y(loc_y)) / (y(loc_y+1) - y(loc_y))             
        
        u1 = val(loc_x,loc_y)   + ratio_x * ( val(loc_x+1,loc_y  ) - val(loc_x,loc_y  ) )  ! change by the first argument, on loc_y
        u2 = val(loc_x,loc_y+1) + ratio_x * ( val(loc_x+1,loc_y+1) - val(loc_x,loc_y+1) )  ! change by the first argument, on loc_y+1
        
        interp = u1 + ratio_y * ( u2 - u1 )

    end subroutine linear_interp2

    ! ####################################################################################################################
    ! 2-dimension linear interpolation, unevenly-spaced grid #############################################################
    ! ####################################################################################################################
    subroutine linear_interp2_uneven(x, xi, y, yi, val, interp)

        implicit none

        double precision, intent(in)    :: x(:)        ! array of grid points
        double precision, intent(in)    :: xi          ! value on the grid line
        double precision, intent(in)    :: y(:)        ! array of grid points
        double precision, intent(in)    :: yi          ! value on the grid line
        double precision, intent(in)    :: val(:,:)    ! array of values corresponding to grid points (x,y)
        double precision, intent(out)   :: interp      ! output, interpolated value
        
        double precision :: u1,u2                      ! intermediate values in interpolation
        
        double precision :: step_x                     ! step size of grid points
        double precision :: gap_x                      ! gap between xi and the lower bound of x
        double precision :: ratio_x                    ! ratio of xi between x(loc) and x(loc+1)
        integer :: loc_x                               ! index for largest grid point among less than x
        
        double precision :: step_y                     ! step size of grid points
        double precision :: gap_y                      ! gap between xi and the lower bound of y
        double precision :: ratio_y                    ! ratio of xi between y(loc) and y(loc+1)
        integer :: loc_y                               ! index for largest grid point among less than y
        
        integer :: i

        integer nx, ny
        nx = size(x)                                   ! number of grid points
        ny = size(y)                                   ! number of grid points

        do i=2, nx-1
            if ( xi <= x(i) ) then
                loc_x = i-1
                exit
            end if
            loc_x = nx-1
        end do
        ratio_x = (xi - x(loc_x)) / (x(loc_x+1) - x(loc_x))

        do i=2, ny-1
            if ( yi <= y(i) ) then
                loc_y = i-1
                exit
            end if
            loc_y = ny-1
        end do
        ratio_y = (yi - y(loc_y)) / (y(loc_y+1) - y(loc_y))
        
        u1 = val(loc_x,loc_y  ) + ratio_x * ( val(loc_x+1,loc_y  ) - val(loc_x,loc_y  ) )  ! change by the first argument, on loc_y
        u2 = val(loc_x,loc_y+1) + ratio_x * ( val(loc_x+1,loc_y+1) - val(loc_x,loc_y+1) )  ! change by the first argument, on loc_y+1

        interp  = u1 + ratio_y * (u2 - u1)

    end subroutine linear_interp2_uneven

    ! ####################################################################################################################
    ! 3-dimension linear interpolation ###################################################################################
    ! ####################################################################################################################
    subroutine linear_interp3(x, xi, y, yi, z, zi, val, interp)   ! Note: works only for evenly-spaced grid!

        implicit none

        double precision, intent(in)    :: x(:)        ! array of grid points
        double precision, intent(in)    :: xi          ! value on the grid line
        double precision, intent(in)    :: y(:)        ! array of grid points
        double precision, intent(in)    :: yi          ! value on the grid line
        double precision, intent(in)    :: z(:)        ! array of grid points
        double precision, intent(in)    :: zi          ! value on the grid line
        double precision, intent(in)    :: val(:,:,:)  ! array of values corresponding to grid points (x,y)
        double precision, intent(out)   :: interp      ! output, interpolated value
        
        double precision :: u11,u21,u12,u22            ! intermediate values in interpolation
        double precision :: w1,w2                      ! intermediate values in interpolation
        
        double precision :: step_x                     ! step size of grid points
        double precision :: gap_x                      ! gap between xi and the lower bound of x
        double precision :: ratio_x                    ! ratio of xi between x(loc) and x(loc+1)
        integer :: loc_x                               ! index for largest grid point among less than x
        
        double precision :: step_y                     ! step size of grid points
        double precision :: gap_y                      ! gap between xi and the lower bound of y
        double precision :: ratio_y                    ! ratio of xi between y(loc) and y(loc+1)
        integer :: loc_y                               ! index for largest grid point among less than y
        
        double precision :: step_z                     ! step size of grid points
        double precision :: gap_z                      ! gap between xi and the lower bound of z
        double precision :: ratio_z                    ! ratio of xi between z(loc) and z(loc+1)
        integer :: loc_z                               ! index for largest grid point among less than z

        integer nx, ny, nz
        nx = size(x)                                   ! number of grid points
        ny = size(y)                                   ! number of grid points
        nz = size(z)                                   ! number of grid points

        step_x  = x(2) - x(1)                                     
        gap_x   = xi - x(1)                                       
        loc_x   = min( nx - 1, max(1, floor(gap_x/step_x) + 1) )    
        ratio_x = (xi - x(loc_x)) / (x(loc_x+1) - x(loc_x))

        step_y  = y(2) - y(1)                                     
        gap_y   = yi - y(1)                                       
        loc_y   = min( ny - 1, max(1, floor(gap_y/step_y) + 1) )    
        ratio_y = (yi - y(loc_y)) / (y(loc_y+1) - y(loc_y))
        
        step_z  = z(2) - z(1)                                     
        gap_z   = zi - z(1)                                       
        loc_z   = min( nz - 1, max(1, floor(gap_z/step_z) + 1) )    
        ratio_z = (zi - z(loc_z)) / (z(loc_z+1) - z(loc_z))

        u11 = val(loc_x,loc_y  ,loc_z) + ratio_x * ( val(loc_x+1,loc_y  ,loc_z) - val(loc_x,loc_y  ,loc_z) )  ! change by the first argument, on loc_y  , loc_z
        u21 = val(loc_x,loc_y+1,loc_z) + ratio_x * ( val(loc_x+1,loc_y+1,loc_z) - val(loc_x,loc_y+1,loc_z) )  ! change by the first argument, on loc_y+1, loc_z
        w1  = u11 + ratio_y * (u21 - u11)

        u12 = val(loc_x,loc_y  ,loc_z+1) + ratio_x * ( val(loc_x+1,loc_y  ,loc_z+1) - val(loc_x,loc_y  ,loc_z+1) )  ! change by the first argument, on loc_y  , loc_z+1
        u22 = val(loc_x,loc_y+1,loc_z+1) + ratio_x * ( val(loc_x+1,loc_y+1,loc_z+1) - val(loc_x,loc_y+1,loc_z+1) )  ! change by the first argument, on loc_y+1, loc_z+1
        w2  = u12 + ratio_y * (u22 - u12)
        
        interp = w1 + ratio_z * ( w2 - w1 )

    end subroutine linear_interp3

    ! ####################################################################################################################
    ! 3-dimension linear interpolation, unevenly-spaced grid #############################################################
    ! ####################################################################################################################
    subroutine linear_interp3_uneven(x, xi, y, yi, z, zi, val, interp)

        implicit none

        double precision, intent(in)    :: x(:)        ! array of grid points
        double precision, intent(in)    :: xi          ! value on the grid line
        double precision, intent(in)    :: y(:)        ! array of grid points
        double precision, intent(in)    :: yi          ! value on the grid line
        double precision, intent(in)    :: z(:)        ! array of grid points
        double precision, intent(in)    :: zi          ! value on the grid line
        double precision, intent(in)    :: val(:,:,:)  ! array of values corresponding to grid points (x,y)
        double precision, intent(out)   :: interp      ! output, interpolated value
        
        double precision :: u11,u21,u12,u22            ! intermediate values in interpolation
        double precision :: w1,w2                      ! intermediate values in interpolation
        
        double precision :: step_x                     ! step size of grid points
        double precision :: gap_x                      ! gap between xi and the lower bound of x
        double precision :: ratio_x                    ! ratio of xi between x(loc) and x(loc+1)
        integer :: loc_x                               ! index for largest grid point among less than x
        
        double precision :: step_y                     ! step size of grid points
        double precision :: gap_y                      ! gap between xi and the lower bound of y
        double precision :: ratio_y                    ! ratio of xi between y(loc) and y(loc+1)
        integer :: loc_y                               ! index for largest grid point among less than y
        
        double precision :: step_z                     ! step size of grid points
        double precision :: gap_z                      ! gap between xi and the lower bound of z
        double precision :: ratio_z                    ! ratio of xi between z(loc) and z(loc+1)
        integer :: loc_z                               ! index for largest grid point among less than z
        integer :: i

        integer nx, ny, nz
        nx = size(x)                                   ! number of grid points
        ny = size(y)                                   ! number of grid points
        nz = size(z)                                   ! number of grid points

        do i=2, nx-1
            if ( xi <= x(i) ) then
                loc_x = i-1
                exit
            end if
            loc_x = nx-1
        end do
        ratio_x = (xi - x(loc_x)) / (x(loc_x+1) - x(loc_x))

        do i=2, ny-1
            if ( yi <= y(i) ) then
                loc_y = i-1
                exit
            end if
            loc_y = ny-1
        end do
        ratio_y = (yi - y(loc_y)) / (y(loc_y+1) - y(loc_y))
        
        do i=2, nz-1
            if ( zi <= z(i) ) then
                loc_z = i-1
                exit
            end if
            loc_z = nz-1
        end do
        ratio_z = (zi - z(loc_z)) / (z(loc_z+1) - z(loc_z))

        u11 = val(loc_x,loc_y  ,loc_z) + ratio_x * ( val(loc_x+1,loc_y  ,loc_z) - val(loc_x,loc_y  ,loc_z) )  ! change by the first argument, on loc_y  , loc_z
        u21 = val(loc_x,loc_y+1,loc_z) + ratio_x * ( val(loc_x+1,loc_y+1,loc_z) - val(loc_x,loc_y+1,loc_z) )  ! change by the first argument, on loc_y+1, loc_z
        w1  = u11 + ratio_y * (u21 - u11)

        u12 = val(loc_x,loc_y  ,loc_z+1) + ratio_x * ( val(loc_x+1,loc_y  ,loc_z+1) - val(loc_x,loc_y  ,loc_z+1) )  ! change by the first argument, on loc_y  , loc_z+1
        u22 = val(loc_x,loc_y+1,loc_z+1) + ratio_x * ( val(loc_x+1,loc_y+1,loc_z+1) - val(loc_x,loc_y+1,loc_z+1) )  ! change by the first argument, on loc_y+1, loc_z+1
        w2  = u12 + ratio_y * (u22 - u12)
        
        interp = w1 + ratio_z * ( w2 - w1 )

    end subroutine linear_interp3_uneven

    ! ####################################################################################################################
    ! 4-dimension linear interpolation ###################################################################################
    ! ####################################################################################################################
    subroutine linear_interp4(x, xi, y, yi, z, zi, w, wi, val, interp)   ! Note: works only for evenly-spaced grid!

        implicit none

        double precision, intent(in)    :: x(:)        ! array of grid points
        double precision, intent(in)    :: xi          ! value on the grid line
        double precision, intent(in)    :: y(:)        ! array of grid points
        double precision, intent(in)    :: yi          ! value on the grid line
        double precision, intent(in)    :: z(:)        ! array of grid points
        double precision, intent(in)    :: zi          ! value on the grid line
        double precision, intent(in)    :: w(:)        ! array of grid points
        double precision, intent(in)    :: wi          ! value on the grid line
        double precision, intent(in)    :: val(:,:,:,:)! array of values corresponding to grid points (x,y)
        double precision, intent(out)   :: interp      ! output, interpolated value
        
        double precision :: u111,u211,u121,u221        ! intermediate values in interpolation
        double precision :: u112,u212,u122,u222        ! intermediate values in interpolation
        double precision :: v11,v21,v12,v22            ! intermediate values in interpolation
        double precision :: interp1,interp2            ! intermediate values in interpolation
        
        double precision :: step_x                     ! step size of grid points
        double precision :: gap_x                      ! gap between xi and the lower bound of x
        double precision :: ratio_x                    ! ratio of xi between x(loc) and x(loc+1)
        integer :: loc_x                               ! index for largest grid point among less than x
        
        double precision :: step_y                     ! step size of grid points
        double precision :: gap_y                      ! gap between xi and the lower bound of y
        double precision :: ratio_y                    ! ratio of xi between y(loc) and y(loc+1)
        integer :: loc_y                               ! index for largest grid point among less than y
        
        double precision :: step_z                     ! step size of grid points
        double precision :: gap_z                      ! gap between xi and the lower bound of z
        double precision :: ratio_z                    ! ratio of xi between z(loc) and z(loc+1)
        integer :: loc_z                               ! index for largest grid point among less than z

        double precision :: step_w                     ! step size of grid points
        double precision :: gap_w                      ! gap between xi and the lower bound of z
        double precision :: ratio_w                    ! ratio of xi between z(loc) and z(loc+1)
        integer :: loc_w                               ! index for largest grid point among less than z

        integer nx, ny, nz, nw
        nx = size(x)                                   ! number of grid points
        ny = size(y)                                   ! number of grid points
        nz = size(z)                                   ! number of grid points
        nw = size(w)                                   ! number of grid points

        step_x  = x(2) - x(1)                                     
        gap_x   = xi - x(1)                                       
        loc_x   = min( nx - 1, max(1, floor(gap_x/step_x) + 1) )    
        ratio_x = (xi - x(loc_x)) / (x(loc_x+1) - x(loc_x))

        step_y  = y(2) - y(1)                                     
        gap_y   = yi - y(1)                                       
        loc_y   = min( ny - 1, max(1, floor(gap_y/step_y) + 1) )    
        ratio_y = (yi - y(loc_y)) / (y(loc_y+1) - y(loc_y))
        
        step_z  = z(2) - z(1)                                     
        gap_z   = zi - z(1)                                       
        loc_z   = min( nz - 1, max(1, floor(gap_z/step_z) + 1) )    
        ratio_z = (zi - z(loc_z)) / (z(loc_z+1) - z(loc_z))

        step_w  = w(2) - w(1)                                     
        gap_w   = wi - w(1)                                       
        loc_w   = min( nw - 1, max(1, floor(gap_w/step_w) + 1) )    
        ratio_w = (wi - w(loc_w)) / (w(loc_w+1) - w(loc_w))

        u111 = val(loc_x,loc_y  ,loc_z,loc_w) + ratio_x * &
               ( val(loc_x+1,loc_y  ,loc_z,loc_w) - val(loc_x,loc_y  ,loc_z,loc_w) )  ! change by the first argument, on loc_y  , loc_z
        u211 = val(loc_x,loc_y+1,loc_z,loc_w) + ratio_x * &
               ( val(loc_x+1,loc_y+1,loc_z,loc_w) - val(loc_x,loc_y+1,loc_z,loc_w) )  ! change by the first argument, on loc_y+1, loc_z
        v11  = u111 + ratio_y * (u211 - u111)

        u121 = val(loc_x,loc_y  ,loc_z+1,loc_w) + ratio_x * &
               ( val(loc_x+1,loc_y  ,loc_z+1,loc_w) - val(loc_x,loc_y  ,loc_z+1,loc_w) )  ! change by the first argument, on loc_y  , loc_z+1
        u221 = val(loc_x,loc_y+1,loc_z+1,loc_w) + ratio_x * &
               ( val(loc_x+1,loc_y+1,loc_z+1,loc_w) - val(loc_x,loc_y+1,loc_z+1,loc_w) )  ! change by the first argument, on loc_y+1, loc_z+1
        v21  = u121 + ratio_y * (u221 - u121)
        
        interp1 = v11 + ratio_z * ( v21 - v11 )

        u112 = val(loc_x,loc_y  ,loc_z,loc_w+1) + ratio_x * &
               ( val(loc_x+1,loc_y  ,loc_z,loc_w+1) - val(loc_x,loc_y  ,loc_z,loc_w+1) )  ! change by the first argument, on loc_y  , loc_z
        u212 = val(loc_x,loc_y+1,loc_z,loc_w+1) + ratio_x * &
               ( val(loc_x+1,loc_y+1,loc_z,loc_w+1) - val(loc_x,loc_y+1,loc_z,loc_w+1) )  ! change by the first argument, on loc_y+1, loc_z
        v12  = u112 + ratio_y * (u212 - u112)

        u122 = val(loc_x,loc_y  ,loc_z+1,loc_w+1) + ratio_x * &
               ( val(loc_x+1,loc_y  ,loc_z+1,loc_w+1) - val(loc_x,loc_y  ,loc_z+1,loc_w+1) )  ! change by the first argument, on loc_y  , loc_z+1
        u222 = val(loc_x,loc_y+1,loc_z+1,loc_w+1) + ratio_x * &
               ( val(loc_x+1,loc_y+1,loc_z+1,loc_w+1) - val(loc_x,loc_y+1,loc_z+1,loc_w+1) )  ! change by the first argument, on loc_y+1, loc_z+1
        v22  = u122 + ratio_y * (u222 - u122)

        interp2 = v12 + ratio_z * ( v22 - v12 )

        interp  = interp1 + ratio_w * ( interp2 - interp1 )

    end subroutine linear_interp4

    ! ####################################################################################################################
    ! 4-dimension linear interpolation, unevenly-spaced grid #############################################################
    ! ####################################################################################################################
    subroutine linear_interp4_uneven(x, xi, y, yi, z, zi, w, wi, val, interp)

        implicit none

        double precision, intent(in)    :: x(:)        ! array of grid points
        double precision, intent(in)    :: xi          ! value on the grid line
        double precision, intent(in)    :: y(:)        ! array of grid points
        double precision, intent(in)    :: yi          ! value on the grid line
        double precision, intent(in)    :: z(:)        ! array of grid points
        double precision, intent(in)    :: zi          ! value on the grid line
        double precision, intent(in)    :: w(:)        ! array of grid points
        double precision, intent(in)    :: wi          ! value on the grid line
        double precision, intent(in)    :: val(:,:,:,:)! array of values corresponding to grid points (x,y)
        double precision, intent(out)   :: interp      ! output, interpolated value
        
        double precision :: u111,u211,u121,u221        ! intermediate values in interpolation
        double precision :: u112,u212,u122,u222        ! intermediate values in interpolation
        double precision :: v11,v21,v12,v22            ! intermediate values in interpolation
        double precision :: interp1,interp2            ! intermediate values in interpolation
        
        double precision :: step_x                     ! step size of grid points
        double precision :: gap_x                      ! gap between xi and the lower bound of x
        double precision :: ratio_x                    ! ratio of xi between x(loc) and x(loc+1)
        integer :: loc_x                               ! index for largest grid point among less than x
        
        double precision :: step_y                     ! step size of grid points
        double precision :: gap_y                      ! gap between xi and the lower bound of y
        double precision :: ratio_y                    ! ratio of xi between y(loc) and y(loc+1)
        integer :: loc_y                               ! index for largest grid point among less than y
        
        double precision :: step_z                     ! step size of grid points
        double precision :: gap_z                      ! gap between xi and the lower bound of z
        double precision :: ratio_z                    ! ratio of xi between z(loc) and z(loc+1)
        integer :: loc_z                               ! index for largest grid point among less than z

        double precision :: step_w                     ! step size of grid points
        double precision :: gap_w                      ! gap between xi and the lower bound of z
        double precision :: ratio_w                    ! ratio of xi between z(loc) and z(loc+1)
        integer :: loc_w                               ! index for largest grid point among less than z
        integer :: i

        integer nx, ny, nz, nw
        nx = size(x)                                   ! number of grid points
        ny = size(y)                                   ! number of grid points
        nz = size(z)                                   ! number of grid points
        nw = size(w)                                   ! number of grid points

        do i=2, nx-1
            if ( xi <= x(i) ) then
                loc_x = i-1
                exit
            end if
            loc_x = nx-1
        end do
        ratio_x = (xi - x(loc_x)) / (x(loc_x+1) - x(loc_x))

        do i=2, ny-1
            if ( yi <= y(i) ) then
                loc_y = i-1
                exit
            end if
            loc_y = ny-1
        end do
        ratio_y = (yi - y(loc_y)) / (y(loc_y+1) - y(loc_y))
        
        do i=2, nz-1
            if ( zi <= z(i) ) then
                loc_z = i-1
                exit
            end if
            loc_z = nz-1
        end do
        ratio_z = (zi - z(loc_z)) / (z(loc_z+1) - z(loc_z))

        do i=2, nw-1
            if ( wi <= w(i) ) then
                loc_w = i-1
                exit
            end if
            loc_w = nw-1
        end do
        ratio_w = (wi - w(loc_w)) / (w(loc_w+1) - w(loc_w))

        u111 = val(loc_x,loc_y  ,loc_z,loc_w) + ratio_x * &
               ( val(loc_x+1,loc_y  ,loc_z,loc_w) - val(loc_x,loc_y  ,loc_z,loc_w) )  ! change by the first argument, on loc_y  , loc_z
        u211 = val(loc_x,loc_y+1,loc_z,loc_w) + ratio_x * &
               ( val(loc_x+1,loc_y+1,loc_z,loc_w) - val(loc_x,loc_y+1,loc_z,loc_w) )  ! change by the first argument, on loc_y+1, loc_z
        v11  = u111 + ratio_y * (u211 - u111)

        u121 = val(loc_x,loc_y  ,loc_z+1,loc_w) + ratio_x * &
               ( val(loc_x+1,loc_y  ,loc_z+1,loc_w) - val(loc_x,loc_y  ,loc_z+1,loc_w) )  ! change by the first argument, on loc_y  , loc_z+1
        u221 = val(loc_x,loc_y+1,loc_z+1,loc_w) + ratio_x * &
               ( val(loc_x+1,loc_y+1,loc_z+1,loc_w) - val(loc_x,loc_y+1,loc_z+1,loc_w) )  ! change by the first argument, on loc_y+1, loc_z+1
        v21  = u121 + ratio_y * (u221 - u121)
        
        interp1 = v11 + ratio_z * ( v21 - v11 )

        u112 = val(loc_x,loc_y  ,loc_z,loc_w+1) + ratio_x * &
               ( val(loc_x+1,loc_y  ,loc_z,loc_w+1) - val(loc_x,loc_y  ,loc_z,loc_w+1) )  ! change by the first argument, on loc_y  , loc_z
        u212 = val(loc_x,loc_y+1,loc_z,loc_w+1) + ratio_x * &
               ( val(loc_x+1,loc_y+1,loc_z,loc_w+1) - val(loc_x,loc_y+1,loc_z,loc_w+1) )  ! change by the first argument, on loc_y+1, loc_z
        v12  = u112 + ratio_y * (u212 - u112)

        u122 = val(loc_x,loc_y  ,loc_z+1,loc_w+1) + ratio_x * &
               ( val(loc_x+1,loc_y  ,loc_z+1,loc_w+1) - val(loc_x,loc_y  ,loc_z+1,loc_w+1) )  ! change by the first argument, on loc_y  , loc_z+1
        u222 = val(loc_x,loc_y+1,loc_z+1,loc_w+1) + ratio_x * &
               ( val(loc_x+1,loc_y+1,loc_z+1,loc_w+1) - val(loc_x,loc_y+1,loc_z+1,loc_w+1) )  ! change by the first argument, on loc_y+1, loc_z+1
        v22  = u122 + ratio_y * (u222 - u122)

        interp2 = v12 + ratio_z * ( v22 - v12 )

        interp  = interp1 + ratio_w * ( interp2 - interp1 )

    end subroutine linear_interp4_uneven

    ! ####################################################################################################################
    ! create a histogram #################################################################################################
    ! ####################################################################################################################
    subroutine histogram(num, max, min, data_vec, bins, count)
    
    ! num  : number of bins
    ! max  : max of the range
    ! min  : min of the range
    ! data_vec : data vector
    ! bins  : num+1 grid for bins, from min to max
    ! count : count in each bin
    
    integer, intent(in) :: num
    double precision, intent(in) :: max, min, data_vec(:)
    double precision :: bin_size

    double precision, intent(out) :: bins(num+1)
    integer, intent(out) :: count(num)
    integer :: i,j
    
    count = 0
    bin_size  = (max - min) / num
    bins(1) = min
    
    do i=1, num
        bins(i+1) = min + bin_size*i
    end do
    
    do i=1, size(data_vec)  ! for each data
    
        do j=2, num+1       ! check which bin it is contained in
            if ( data_vec(i) < bins(j) ) then
               count(j-1) = count(j-1) + 1
               exit
            end if
        end do
        
    end do
    
    end subroutine histogram

    ! ####################################################################################################################
    ! binomial distribution ##############################################################################################
    ! ####################################################################################################################
    function binom(k,n,p)
    
    integer :: i,n,k,ncr
    double precision :: p,binom,log_ncr
    
    log_ncr = 0d0
    
    if ( n < 0 ) then
        write(*,*) n, "n is negative!"
    end if
    if ( k < 0 ) then
        write(*,*) k, "k is negative!"
    end if
    if ( n < k ) then
        write(*,*) "n < k !"
    end if
        
    if ( n == 0 ) then
    
        binom = 1d0
        
    else
            
        if ( k == 0 ) then
            ncr = 1
        elseif ( n-k == 0 ) then
            ncr = 1
        else
            do i = n-k+1, n
                log_ncr = log_ncr + dlog10(dble(i))
            end do
            do i = 1, k
                log_ncr = log_ncr - dlog10(dble(i))
            end do
            ncr = nint(10d0**log_ncr)
            
        end if

        binom = ncr * p**k * (1d0-p)**(n-k)
    
    end if
    
    end function binom

    ! ####################################################################################################################
    ! linear numerical derivatives #######################################################################################
    ! ####################################################################################################################
    subroutine num_drv(x, xi, val, drv)

        implicit none

        double precision, intent(in)    :: x(:)        ! array of grid points
        double precision, intent(in)    :: xi          ! value on the grid line
        double precision, intent(in)    :: val(:)      ! array of values corresponding to grid points x
        double precision, intent(out)   :: drv         ! output, numerical derivative
                
        integer :: loc_x                               ! index for the largest grid point lower than x
        integer :: i
        integer nx
        
        nx = size(x)                                   ! number of grid points

        do i=2, nx-1
            if ( xi <= x(i) ) then
                loc_x = i-1
                exit
            end if
            loc_x = nx-1
        end do
        
        drv = (val(loc_x+1) - val(loc_x)) / (x(loc_x+1) - x(loc_x))

    end subroutine num_drv
    
    ! ###########################################################################################################
    ! trapezoidal integration ###################################################################################
    ! ###########################################################################################################
    subroutine trapz_integ(x,y,value)
    
    ! input x: grid points across the range of integration
    !          need to be strictly increasing; no need to be evenly-spaced
    ! input y: values corresponding to each grid on the range

        implicit none

        double precision, intent(in)  :: x(:),y(:)
        double precision, intent(out) :: value
        double precision :: trapz
        integer :: num_x,num_y,i
        
        num_x = size(x)
        num_y = size(y)
        if ( num_x /= num_y ) then
            write(*,*) "size of x and y do not match!"
        end if

        value = 0d0
        do i=2, num_x-1
            trapz = ( y(i+1)+y(i) ) * ( x(i+1) - x(i) ) / 2d0
            value = value + trapz
        end do

    end subroutine trapz_integ

    ! ###########################################################################################################
    ! resize 2-dimensional matrix with real values ##############################################################
    ! ###########################################################################################################
    subroutine resize_2_real(matrix,x,y)
    
    ! When the original matrix is truncated, the contents are preserved.
    ! When the original matrix is augmented, the contents are lost and set to zeros.

    implicit none

        double precision, allocatable :: matrix(:,:)
        integer, intent(in) :: x,y
        
        deallocate(matrix)
        allocate(matrix(x,y))
        
    end subroutine resize_2_real

    ! ####################################################################################################################
    ! normal cumulative function from asa066 #############################################################################
    ! ####################################################################################################################

subroutine normp ( z, p, q, pdf )

!*****************************************************************************80
!
!! NORMP computes the cumulative density of the standard normal distribution.
!
!  Discussion:
!
!    This is algorithm 5666 from Hart, et al.
!
!  Modified:
!
!    13 January 2008
!
!  Author:
!
!    Original FORTRAN77 version by Alan Miller.
!    FORTRAN90 version by John Burkardt.
!
!  Reference:
!
!    John Hart, Ward Cheney, Charles Lawson, Hans Maehly,
!    Charles Mesztenyi, John Rice, Henry Thacher,
!    Christoph Witzgall,
!    Computer Approximations,
!    Wiley, 1968,
!    LC: QA297.C64.
!
!  Parameters:
!
!    Input, real ( kind = 8 ) Z, divides the real line into two
!    semi-infinite intervals, over each of which the standard normal
!    distribution is to be integrated.
!
!    Output, real ( kind = 8 ) P, Q, the integrals of the standard normal
!    distribution over the intervals ( - Infinity, Z] and
!    [Z, + Infinity ), respectively.
!
!    Output, real ( kind = 8 ) PDF, the value of the standard normal
!    distribution at Z.
!
  implicit none

  real ( kind = 8 ) :: cutoff = 7.071D+00
  real ( kind = 8 ) expntl
  real ( kind = 8 ) p
  real ( kind = 8 ) :: p0 = 220.2068679123761D+00
  real ( kind = 8 ) :: p1 = 221.2135961699311D+00
  real ( kind = 8 ) :: p2 = 112.0792914978709D+00
  real ( kind = 8 ) :: p3 = 33.91286607838300D+00
  real ( kind = 8 ) :: p4 = 6.373962203531650D+00
  real ( kind = 8 ) :: p5 = 0.7003830644436881D+00
  real ( kind = 8 ) :: p6 = 0.03526249659989109D+00
  real ( kind = 8 ) pdf
  real ( kind = 8 ) q
  real ( kind = 8 ) :: q0 = 440.4137358247522D+00
  real ( kind = 8 ) :: q1 = 793.8265125199484D+00
  real ( kind = 8 ) :: q2 = 637.3336333788311D+00
  real ( kind = 8 ) :: q3 = 296.5642487796737D+00
  real ( kind = 8 ) :: q4 = 86.78073220294608D+00
  real ( kind = 8 ) :: q5 = 16.06417757920695D+00
  real ( kind = 8 ) :: q6 = 1.755667163182642D+00
  real ( kind = 8 ) :: q7 = 0.08838834764831844D+00
  real ( kind = 8 ) :: root2pi = 2.506628274631001D+00
  real ( kind = 8 ) z
  real ( kind = 8 ) zabs

  zabs = abs ( z )
!
!  37 < |Z|.
!
  if ( 37.0D+00 < zabs ) then

    pdf = 0.0D+00
    p = 0.0D+00
!
!  |Z| <= 37.
!
  else

    expntl = exp ( - 0.5D+00 * zabs * zabs )
    pdf = expntl / root2pi
!
!  |Z| < CUTOFF = 10 / sqrt(2).
!
    if ( zabs < cutoff ) then

      p = expntl * (((((( &
          p6   * zabs &
        + p5 ) * zabs &
        + p4 ) * zabs &
        + p3 ) * zabs &
        + p2 ) * zabs &
        + p1 ) * zabs &
        + p0 ) / ((((((( &
          q7   * zabs &
        + q6 ) * zabs &
        + q5 ) * zabs &
        + q4 ) * zabs &
        + q3 ) * zabs &
        + q2 ) * zabs &
        + q1 ) * zabs &
      + q0 )
!
!  CUTOFF <= |Z|.
!
    else

      p = pdf / ( &
        zabs + 1.0D+00 / ( &
        zabs + 2.0D+00 / ( &
        zabs + 3.0D+00 / ( &
        zabs + 4.0D+00 / ( &
        zabs + 0.65D+00 )))))

    end if

  end if

  if ( z < 0.0D+00 ) then
    q = 1.0D+00 - p
  else
    q = p
    p = 1.0D+00 - q
  end if

  return
end 

    ! ####################################################################################################################
    ! function minimization from nms #####################################################################################
    ! ####################################################################################################################
    function fmin ( a, b, f, tol )

    !*****************************************************************************80
    !
    !  FMIN seeks a minimizer of a scalar function of a scalar variable.
    !
    !  Discussion:
    !
    !    FMIN seeks an approximation to the point where F attains a minimum
    !    value (strictly) inside the interval (A,B).
    !
    !    The method used is a combination of golden section search and
    !    successive parabolic interpolation.  Convergence is never much
    !    slower than that for a Fibonacci search.  If F has a continuous
    !    second derivative which is positive at the minimum (which is not
    !    at A or B), then convergence is superlinear, and usually of the
    !    order of about 1.324....
    !
    !    The function F is never evaluated at two points closer together
    !    than EPS * ABS ( XSTAR ) + (TOL/3), where EPS is approximately the
    !    square root of the relative machine precision.  If F is a unimodal
    !    function and the computed values of F are always unimodal when
    !    separated by at least EPS * ABS ( XSTAR ) + (TOL/3), then FMIN
    !    approximates the abcissa of the global minimum of F on the
    !    interval [A, B] with an error less than 3 * EPS * ABS ( FMIN ) + TOL.
    !    If F is not unimodal, then FMIN may approximate a local, but
    !    perhaps non-global, minimum to the same accuracy.
    !
    !    It is worth stating explicitly that this routine will NOT be
    !    able to detect a minimizer that occurs at either initial endpoint
    !    A or B.  If this is a concern to the user, then the user must
    !    either ensure that the initial interval is larger, or to check
    !    the function value at the returned minimizer against the values
    !    at either endpoint.
    !
    !  Modified:
    !
    !    16 October 2004
    !
    !  Reference:
    !
    !    Richard Brent,
    !    Algorithms for Minimization Without Derivatives,
    !    Dover, 2002,
    !    ISBN: 0-486-41998-3,
    !    LC: QA402.5.B74.
    !
    !    David Kahaner, Cleve Moler, Steven Nash,
    !    Numerical Methods and Software,
    !    Prentice Hall, 1989,
    !    ISBN: 0-13-627258-4,
    !    LC: TA345.K34.
    !
    !  Parameters
    !
    !    Input/output, real ( kind = 8 ) A, B.  On input, the left and right
    !    endpoints of the initial interval.  On output, the lower and upper
    !    bounds for the minimizer.  It is required that A < B.
    !
    !    Input, external F, a real function of the form
    !      function f ( x )
    !      real f
    !      real x
    !    which evaluates F(X) for any X in the interval [A,B].
    !
    !    Input, real ( kind = 8 ) TOL, the desired length of the interval of
    !    uncertainty of the final result.  TOL must not be negative.
    !
    !    Output, real ( kind = 8 ) FMIN, the approximate minimizer of the function.
    !
    implicit none

    real ( kind = 8 ) a
    real ( kind = 8 ) b
    real ( kind = 8 ) c
    real ( kind = 8 ) d
    real ( kind = 8 ) e
    real ( kind = 8 ) eps
    real ( kind = 8 ), external :: f
    real ( kind = 8 ) fmin
    real ( kind = 8 ) fu
    real ( kind = 8 ) fv
    real ( kind = 8 ) fw
    real ( kind = 8 ) fx
    real ( kind = 8 ) midpoint
    real ( kind = 8 ) p
    real ( kind = 8 ) q
    real ( kind = 8 ) r
    real ( kind = 8 ) tol
    real ( kind = 8 ) tol1
    real ( kind = 8 ) tol2
    real ( kind = 8 ) u
    real ( kind = 8 ) v
    real ( kind = 8 ) w
    real ( kind = 8 ) x

    if ( b <= a ) then
        write ( *, '(a)' ) ' '
        write ( *, '(a)' ) 'FMIN - Fatal error!'
        write ( *, '(a)' ) '  A < B is required, but'
        write ( *, '(a,g14.6)' ) '  A = ', a
        write ( *, '(a,g14.6)' ) '  B = ', b
        stop
    end if

    c = 0.5D+00 * ( 3.0D+00 - sqrt ( 5.0D+00 ) )
    !
    !  C is the squared inverse of the golden ratio.
    !
    !  EPS is the square root of the relative machine precision.
    !
    eps = sqrt ( epsilon ( eps ) )
    !
    !  Initialization.
    !
    v = a + c * ( b - a )
    w = v
    x = v
    e = 0.0D+00
    fx = f(x)
    fv = fx
    fw = fx
    !
    !  The main loop starts here.
    !
    do

        midpoint = 0.5D+00 * ( a + b )
        tol1 = eps * abs ( x ) + tol / 3.0D+00
        tol2 = 2.0D+00 * tol1
    !
    !  Check the stopping criterion.
    !
        if ( abs ( x - midpoint ) <= ( tol2 - 0.5D+00 * ( b - a ) ) ) then
            exit
        end if
    !
    !  Is golden-section necessary?
    !
        if ( abs ( e ) <= tol1 ) then
            if ( midpoint <= x ) then
                e = a - x
            else
                e = b - x
            end if

        d = c * e
    !
    !  Consider fitting a parabola.
    !
        else

            r = ( x - w ) * ( fx - fv )
            q = ( x - v ) * ( fx - fw )
            p = ( x - v ) * q - ( x - w ) * r
            q = 2.0D+00 * ( q - r )
            if ( 0.0D+00 < q ) then
                p = -p
            end if
            q = abs ( q )
            r = e
            e = d
    !
    !  Choose a golden-section step if the parabola is not advised.
    !
            if ( &
               ( abs ( 0.5D+00 * q * r ) <= abs ( p ) ) .or. &
               ( p <= q * ( a - x ) ) .or. &
               ( q * ( b - x ) <= p ) ) then

                if ( midpoint <= x ) then
                    e = a - x
                else
                    e = b - x
                end if

            d = c * e
    !
    !  Choose a parabolic interpolation step.
    !
            else

                d = p / q
                u = x + d

                if ( ( u - a ) < tol2 ) then
                    d = sign ( tol1, midpoint - x )
                end if

                if ( ( b - u ) < tol2 ) then
                    d = sign ( tol1, midpoint - x )
                end if

            end if

        end if
    !
    !  F must not be evaluated too close to X.
    !
        if ( tol1 <= abs ( d ) ) then
            u = x + d
        end if

        if ( abs ( d ) < tol1 ) then
            u = x + sign ( tol1, d )
        end if

        fu = f(u)
    !
    !  Update the data.
    !
        if ( fu <= fx ) then

            if ( x <= u ) then
                a = x
            else
                b = x
            end if

            v = w
            fv = fw
            w = x
            fw = fx
            x = u
            fx = fu
            cycle

        end if

        if ( u < x ) then
            a = u
        else
            b = u
        end if

        if ( fu <= fw .or. w == x ) then
            v = w
            fv = fw
            w = u
            fw = fu
        else if ( fu <= fv .or. v == x .or. v == w ) then
            v = u
            fv = fu
        end if

    end do

    fmin = x

    return
    
    end function

end module